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The cumulant expansion is used to estimate generalized Lyapunov exponents of the random- 
frequency harmonic oscillator. Three stochastic processes are considered: Gaussian white noise, 
Ornstein-Uhlenbeck, and Poisson shot noise. In some cases, nontrivial numerical difficulties arise. 
These are mostly solved by implementing an appropriate importance-sampling Montecarlo scheme. 
We analyze the relation between random-frequency oscillators and many-particle systems with pair- 
wise interactions like the Lennard- Jones gas. 
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I. INTRODUCTION 

Lyapunov exponents quantify sensitivity to initial con- 
ditions in dynamical systems. The existence of a posi- 
tive Lyapunov exponent implies that trajectories initially 
close in phase space will typically diverge exponentially 
fast in time. In practice, this sets a limit for predict- 
ing the future behavior of the system, because small im- 
precisions in the knowledge of the initial state will be 
amplified at a rate given by the largest Lyapunov expo- 
nent. Even if determinism subsists on a short time scale, 
on longer time windows the system exhibits features of 
randomness. This property lies at the basis of the statis- 
tical description of many-particle deterministic systems. 
Hence the interest in analytical estimates of Lyapunov 
exponents in simple statistical-mechanics models. 

The theory of Lyapunov exponents of hard-ball sys- 
tems has a long history. It started with the pioneering 
work of Krylov [l|, |3| , was rigorously developed by Sinai 
3 and collaborators, and completed (to some extent) by 
van Beijeren, Dorfman and co-workers . The analyt- 
ical calculation of, e.g., the largest Lyapunov exponent 
of a dilute rigid-sphere gas, is based on the fact that the 
dynamics consists of free rectilinear motions interrupted 
by instantaneous elastic collisions 6] ; the expressions so- 
obtained agree quantitatively with the numerical exper- 
iments 

The case of a dilute gas with finite-range interactions 
can be handled in close analogy with the rigid-sphere 
problem: though the collisions are not trivial any more, 
the dynamics is still ruled by occasional pairwise encoun- 
ters [a, • However, when one considers long-range 
interactions (or short-range interactions and high den- 
sities), the theoretical approach must be substantially 
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modified. 

In the general case we must deal with the full system 
of coupled differential equations that govern the evolu- 
tion of multidimensional tangent vectors ?y(t). Consider 
for instance a gas of N particles in three dimensions de- 
scribed by the Hamiltonian 
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where q,; and pi, are conjugate position-momentum co- 
ordinates. Assuming m = 1, tangent vectors evolve ac- 
cording to 



r 



(2) 



(dot meaning time derivative), where V is the Hessian 
matrix of the potential V, namely 



dqtdq-j 



(3) 



The Hessian depends explicitly on time because it is cal- 
culated along a reference trajectory q{t). Once initial 
conditions zq = {qo,Po) and rjo have been specified, one 
can find r]{t) from Eq. ([2|). Then the Lyapunov exponent 
A is obtained by calculating the limit |13j | 



where 



A = lim A (<; ^0,770) 

t^OO 



A (i; 2:0,770) = -ln|?7(t;zQ,?7o) 



(4) 



(5) 



The finite-time Lyapunov exponent A(t;zoj'7o) depends 
on the initial conditions zq and 770. However, assuming 
ergodicity on the energy-shell, A becomes independent of 
Zq, which can then be chosen randomly, e.g., according 
to the microcanonical distribution. There will also be 
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no dependence on initial tangent vectors, because if 770 is 
also chosen randomly, it will always have a non-zero com- 
ponent along the most expanding direction. In spite of 
being redundant, the averaging over zq and rjQ permits to 
treat equations ([2]) foi'mally as a system of stochastic dif- 
ferential equations (l4|. So, in this "stochastic" approach 
one attempts the analytical estimation of the average 



A = lim i(ln|?7(t;zo,'7o)|) 

t— >oo t 
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This is a hard task, however. It is much simpler to eval- 

(7) 



simp ler 1 

uate the generalized Lyapunov exponent [l5l [l6j 
1 



A2 = lim ^ln(|?7(i;2;o,?7o)n 



and assume it approximately coincides with the standard 
Lyapunov exponent, which is justified in the absence of 
intermittency [l6l |. 

Moreover, if the Hamiltonian of the system can be de- 
composed as some "free" part plus weak interactions, 
then perturbative techniques, like the cumulant expan- 
sion [3 [13, can be invoked. This is essentially the 
approach followed by Barnett et al. fiol- 2TI. Pettini et 
al. dHH, and the present authors [lllig. Though 
there are some differences among the formulations of the 
three groups above, it may be said that the main the- 
oretical conclusion extracted from that body of work is 
the following. As far as A2 is concerned, if one combines 
the cumulant expansion with some kind of isotropy ap- 
proximation (which may be fully justified in some cases), 
the original problem of 6iV differential equations can be 
reduced to a system of only two equations for a "repre- 
sentative" single degree of freedom: 
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In this kind of mean-field approximation, the "curvature" 
K(t) is a scalar stochastic process, whose cumulants can 
be related to the operator cumulants of the Hessian V(t) 
(see, e.g., (H). 

The comparison of theoretical results obtained with 
the cumulant expansion -truncated at the second order- 
versus numerical simulations has met mixed success. The 
agreement is very good for some many-particle systems 
with bounded weak interactions [2^, [23] and for the 
Fermi-Pasta-Ulam system 2J]. On the other side, the 
results for the Id-XY model [l^ and for a dense one- 
component plasma [l^ [2^ are not so satisfactory. 

Anyway, the mentioned tests, which compare theoreti- 
cal estimates for A2 against numerical calculations for A, 
should be taken with some reservations: (a) Pettini et 
al. did not check if the approximate equality A « A2 in- 
deed holds [23 - [23 |. Moreover, their theory includes a fit- 
ting parameter [the correlation time of the process K(t)]. 
Then, it may happen that the theory really agrees with 
the simulations, or, alternatively, it may be the case of 
a disagreement that is compensated by a suitable choice 



of the correlation time, (b) The authors of Refs. |25l - l27j | 
derived Eq. ([8|) from first principles (no fitting parame- 
ters) and verified numerically that A ~ A2 holds in their 
tests. However, they used a simple ("brute force" |2j|) 
Montecarlo s amp ling for doing the average ([7]). And it is 
known (e.g., [23]) that simple samplings tend to produce 
wrong estimates of generalized Lyapunov exponents 



Aq = lim -i-ln(]?7(i;zo,?7o)l'^ 
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The larger the value of Q, the stronger this spurious ef- 
fect. (Consistently, there are no difficulties in the numer- 
ical calculation of the standard A, given that Ag A for 
Q^O.) 

In conclusion: if one wants to assess the quality of the- 
oretical predictions unabiguously, then it is necessary to 
develop trustable Montecarlo algorithms for the calcula- 
tion of Aq. We are not aware of the existence of such 
methods for Hamiltonian many-particle systems. On 
the other side, an importance-sampling ^30i] algorithm 
was recently proposed by Vanneste for calculating Aq in 
stochastic dynamical systems. The algorithm was shown 
to perform efficiently for white noise and Q not too large 

Hi. 

The present work is part of a larger project that aims 
at defining the limits of validity of the cumulant approach 
for the Lyapunov exponent of many-particle Hamiltonian 
systems. We start our investigations with the simpli- 
fied mean-field setting ([5]). This is the simplest possible 
case having the same formal structure as the many-body 
problem. By choosing K,{t) to be a stochastic process we 
shall be able to use importance-sampling in the numeri- 
cal calculations. For several choices of hi{t) we shall both 
analyze the performance of the cumulant expansion and 
test the numerical algorithms. 

It has been argued [2^] that, for typical chaotic many- 
body systems, K,{t) should be close to Gaussian white 
noise; this is the first case we shall consider. For Gaussian 
white-noise the second-order cumulant expansion for A2 
is exact, thus this case is ideally suited for analyzing the 
difficulties that appear in the numerical calculation of A2 
(Sect.HYl. 

Next, we keep the Gaussian and Markov properties but 
allow for finite correlation times, leading to the Ornstein- 
Uhlcnbeck process. In this case we calculate the fourth 
cumulant contribution to A2. This test will give us 
some idea of (i) the convergence rate of the cumulant 
expansion, and (ii) the performance of the importance- 
sampling method for colored noise fSect. |V|) . 

Last we study the situation of K.{t) being Poisson white 
shot-noise. This appears to be the appropriate choice 
for modeling the tangent-vector dynamics in dilute gases 
with short-range interactions. Like in the case of Gaus- 
sian white noise, here we have analytical expressions for 
the generalized exponents A2, A4, Ae, etc. So, this case 
will provide an opportunity for further testing of the nu- 
merical algorithm. At the same time it will be helpful for 
characterizing the distribution of finite-time Lyapunov 
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exponents, e.g., when is this distribution approximately 
Gaussian? (Sect. EH). 

Section contains a short review of the cumulant ex- 
pansion as appUed to the determination of some gen- 
erahzed Lyapunov exponents. In Sect. IIIII we describe 
the three Montecarlo methods considered in this paper: 
simple, simple-Gaussian, and importance-sampling. Sec- 
tion |VlT] presents a summary of our results and the final 
remarks. 

Before proceeding to the bulk of the paper let us com- 
ment that the random oscillator of Eq. is formally 
equivalent to the Schrodinger equation for a particle in 
a disordered potential (Anderson localization problem in 
one dimension) . Thus many useful results related to ran- 
dom oscillators can be found in the condensed-matter 
literature jsMil. 



Ellipsis stand for third and higher cumulants (some ex- 
plicit expressions can be found in [18]). The exponent A2 
is related to the eigenvalue of K that has the largest real 
part: 

\2 ^ ^iTiax^ {ki,k2,k3} , (15) 

with ki the eigenvalues of K. 

Starting from the evolution equations for higher order 
products [analogous to ([TT|) ] and repeating the same steps 
above, one can derive the corresponding expressions for 
A4, Ae, etc. Of course, the algebraic difficulties increase 
with the order of the exponent. 



III. NUMERICAL METHODS 



II. CUMULANT EXPANSION FOR THE KUBO 
OSCILLATOR 

Equation ([8]) describes a harmonic oscillator with a 
random frequency ui such that iu'^ — k (Kubo oscillator) . 
It is worth extending this model a bit to account for the 
possibility of damping, i.e., we shall consider an oscillator 
described by the first-order equations 

q = P, 

p + ap + Kq = 0. (10) 

Let us make the identifications q ~ rji, p — ri2 [59] . Then, 
putting Of = we recover ([8]). 

Some analytical results for the Lyapunov exponent of 
the Kubo oscillator pU|) can be found in the literature 
(see, e.g., [iol - li^ ). Here we shall concentrate on the gen- 
eralized exponent A2. For this purpose we must consider 
the dynamics of second order products: 
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Let us assume that both parameters a and k are station- 
ary stochastic processes. If fluctuations are small enough 
(in a sense that will be discussed later), one can obtain 
dynamical equations for the second-order averages using 
the cumulant expansion [l^ . Splitting the stochastic ma- 
trix as an average plus fluctuations: 

B(t) = Bo + Bi(i), (12) 

it can be shown that for long times one has [l^ : 




(13) 

where K is the 3x3 matrix given by the operator cumu- 
lant expansion 14] 

/>oo 

K = Bo+ / (Bi(T)e^«"Bi(0))e-^«"dT+... . (14) 
Jo 



The numerical evolution of Eqs. (jlO[) was performed by 
means of the Euler algorithm with time step dt = 10~^ 
(some higher-order algorithms |43|] were tested, but did 
not lead to substantial improvements). A set of trajec- 
tories is generated by randomly choosing (go,Po); (^{t) 
and K{t). For each trajectory we computed the norm 
\ri{t)\ = ^ q^ + as a function of time. The Lyapunov 
exponent is then approximated by the average of finite- 
time exponents: 



A«(A(t;Co))-(7ln|»7(<;Co)|) 



(16) 



The operation (• ■ ■ ) means averaging over a certain num- 
ber of realizations of the pseudorandom variables (com- 
pactly denoted by Co) that determine the trajectories. 
Time t must be large enough to assure the convergence 
of the average to the desired precision. 

In principle we could use the same scheme as before 
for estimating generalized exponents, i.e.. 



^MI.(t;Co)|«)^^ln(e 



QtA(t:Co) 



(17) 



the last equality following from (1161) . However such a sim- 
ple averaging tends to undestimate rare events. Hence, 
spurious results are expected whenever the distribution 
P(At) does not decay fast enough [2^, [H, H^] ■ A some- 
what better alternative is, instead of straightforward av- 
eraging, to estimate the generalized exponent from the 
first terms of the series 



(Qt)"-1 



(18) 



where the nth-order cumulants of P {\t) [43 • In 

principle, these cumulants could be estimated numeri- 
cally. However, for the samples we considered, third and 
higher cumulants are typically rather unstable [4^ . So, 
it is practically impossible to assess the convergence of 
the expansion (|18p . For this reason, cumulants k„, with 



4 



n > 3, will not be included in our calculations. Thus we 
arrive at 

^Q~X+^QtK2{t). (19) 

[If P{Xt) is Gaussian, this expression is exact.] We call 
the procedure leading to Eq. (|19|) simple Gaussian aver- 
aging. From Eq. (|19p one can derive approximate expres- 
sions for the standard Lyapunov exponent, the simplest 
one being 

A « 2A2 - A4 . (20) 

Conversely, when A, A2, A4 are known, the deviation from 
equality in the formula above provides a measure of the 
non-Gaussianity of P{Xt). 

When the tail of P{Xt) is essential for the determina- 
tion of Aq and it is not Gaussian, the approximations (jl7p 
and (IT9I) are bound to fail. In this case one must resort 
to numerical methods capable of sampling the relevant 
part of the distribution P(Ai). The importance-sampling 
Monte Carlo algorithm recently proposed by Vanneste 
is especially suited for our needs. The algorithm, both 
efficient and easy to implement, uses a simple random 
resampling step: those trajectories which contribute the 
most (least) to the average are cloned (pruned) with a 
large probability [29j . 

Having presented the theory and the numerical meth- 
ods, we are ready to proceed with the comparisons. 



After substitution into Eq. ([T4)) we readily obtain 
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The generalized exponent A2 can now be calculated from 
Eq. ([T5l) . A closed expression for the standard L yap unov 
exponent was derived by Mallick and Peyneau [405. As 
an example, we display in Fig. [T] analytical and numeri- 
cal results for both exponents. Given that the theoretical 
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IV. GAUSSIAN WHITE NOISE 



Only when the matrix stochastic process Bi is Gaus- 
sian and delta-correlated the cumulant expansion does 
stop at the second order, i.e., Eq. (fT4)) without the el- 
lipsis becomes exact 18|. This is the case we consider 
now. 

(Stochastic differential equations with multiplicative 
white noise will always be taken in the Stratonovich 
sense.) 



FIG. 1: Lyapunov exponents versus noise strength for the 
harmonic oscillator with random frequency. Symbols cor- 
respond to numerical results for A (small circles) and A2 
(large circles: importance-sampling; triangles: simple Gaus- 
sian sampling). Two values of the damping constant were 
used: a = (hollow symbols) and a = 1 (filled symbols). In 
both cases kq = 1. In all cases we averaged over 10"^ trajecto- 
ries. Resampling time was set to t^cs = 1-0. Lines correspond 
to exact theoretical expressions. 



A. Random frequency 

Let us first study the situation where the damping a 
is a constant and 



(21) 



where ^ (t) is a zero- mean Gaussian white noise with cor- 
relation function 



With these definitions one has 
2 



(22) 
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expressions are exact, this comparison constitutes a rigor- 
ous test for the numerical methods. We see that, even for 
relatively small samples, the importance-sampling calcu- 
lation agrees perfectly with the theory. The Gaussian 
sampling, though not perfect, provides a reasonably good 
approximation. 

Clearly both exponents, A2 and A, do not coincide. 
This is to be expected whenever fluctuations in the fre- 
quency/damping are large as compared to their average 
values [3a, 14^ . 

The higher order exponents A2J, with J = 2, 3, . . . are 
obtained by diagonalizing matrices of size 2 J — 1. Such 
matrices describe the evolution of the moments (<;"p'"), 
with m-\-n = 2 J, and have simple analytical expressions 
[35I . IZsj . An example involving a higher order exponent 
will be shown in Sec. I VII 



5 



B. Random damping 

Now we consider a harmonic oscillator with constant 
frequency but in an environment with fluctuating damp- 
ing coefficient 



a{t) = ao + ^(i) , 



(25) 



where ^(t) is again zero-mean Gaussian white noise, with 
correlation given by Eq. 
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FIG. 2: Lyapunov exponents vs. noise strength for the har- 
monic oscillator with random damping. Lines represent exact 
theoretical results. Symbols correspond to numerical calcula- 
tions for A2 (large circles) and A (small circles). Importance- 
sampling Monte Carlo was used in the case of A2. We chose 
two values for the average damping coefficient: ao — 1 (top 
panel), qq = 4 (bottom panel). In both cases k = 1. In all 
cases we averaged over 10^ trajectories. Resampling time was 

set to trcs = 0.2. 

Now the matrix B is decomposed as 



B 



2 

-2ao -2k 

— K 1 — Q!o 
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• 










(26) 



Hence, substitution into ([H)) yields 



K 





-2ao + 2A 
1 




(27) 



Upon diagonalizing K we obtain A2 . Figure [2] presents 
the comparison of numerical and analytical results for A 
and A2 as the noise intensity is varied (exact theoretical 
results for A were extracted from Ref. 41]). 

Concerning the numerical calculation of A2, besides 
noting the excellent agreement with the theory, it must 
be said that the importance-sampling method behaved in 
a very robust way both for random frequency and ran- 
dom damping. Changing sample size, simulation time, 
and resampling time t^es j29i] within reasonable limits did 
not appreciably affect the result for A2. However, if one 
increases ires beyond certain bounds, then the method 
becomes inefficient, as very large samples are necessary 
to guarantee convergence to the correct results. 



V. CORRELATED NOISE 

For white noise fluctuations, either in the frequency or 
in the damping, we have verified in the previous section 
that the theory for A2 is in agreement with numerical 
results, provided the latter are obtained using importance 
sampling. 

Now we shall analyze the effect of introducing noise 
correlations. We consider the case of a random frequency, 
as in Eq. (j21[) . but now the noise is a zero-mean Ornstein- 
Ulhenbeck process, i.e., with correlation function 

mm) = ^ cxp(-|t - t'lM EE a' cxp(-|< - t'\/r) . 

(28) 

For simplicity we set a — and kq = 0. By inserting 
p3| into (|14p . the second-cumulant matrix becomes 



K(2) ^ 




2 


-2A' 



(29) 



Notice that in the limit t — > the white-noise case is 
recovered. 

In the presence of correlations the second-order trun- 
cation of the cumulant expansion (jl4p is not exact. In 
order to improve the theory one must calculate higher cu- 
mulants. For the present case the third cumulant is null. 
Explicit expressions for the fourth cumulant were given 
by Fox [31 , Breuer et al. [1^ , and Tessieri [34| . A some- 
what lengthy calculation (sketched in Appendix |X| leads 
to the following result for the fourth order approximation 
to K: 



K(4) = K( 




(30) 



The comparison between numerical and theoretical re- 
sults for A2 is presented in Fig. [31 For completeness we 
also show numerical calculations of the standard Lya- 
punov exponent, together with an approximate theoret- 
ical expression obtained along the lines of Ref. [1^ (see 
Appendix |B| . 
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FIG. 3: (Color online.) Harmonic oscillator with correlated 
random frequency. Symbols indicate numerical results for A2 
(full circles) and A (hollow circles) as a function of the noise 
amplitude A (averages over 10* trajectories). Parameters are 
a — 0, Ko = and t — 1. Resampling time was set to trcs ~ 
20. Solid lines correspond to theoretical results. For A2 we 
used the cumulant expansion (blue: truncation at the second 
cumulant; dark blue: including the fourth cumulant). An 
approximate analytical expression for A (red) is also shown. 



We see in Fig. [3] that the inclusion of the fourth cu- 
mulant contribution noticeably extends the domain of 
validity of the theory into the region of larger noise am- 
plitudes (with respect to the second-order approxima- 
tion). Higher cumulants can also be calculated, but the 
required effort quickly becomes unbearable. For instance, 
the sixth cumulant demands the calculation of more than 
100 terms (see Appendix E]) . Anyway, the theory being 
perturbative, by increasing the amplitude of the noise 
and/or the correlation time, one eventually arrives at 
a point were the cumulant expansion completely breaks 
down. 

The perturbation parameter controlling the conver- 
gence of the cumulant expansion is the so-called Kubo 
number e. General considerations led van Kampen [l4] to 
conclude that the Kubo number is the product of the am- 
plitude of the fluctuations and the correlation time, that 
is (TT. However, in the present case it is clear that such 
a combination is not adimensional. The correct Kubo 
number is instead 



(31) 



This can be checked explicitly from the second and fourth 
cumulants above. Consider, for instance, the element 
K21, which dominates the Lyapunov exponent for small 
correlation times: 



A- 



13 



AM 



13 



-Ar' 



In the white- noise limit, i.e., r — with A 
Kubo number tends to zero -as it should be. 



fixed. 



(32) 
the 



On the numerical side, we comment that for large noise 
amplitudes the convergence to the limiting values is much 
slower than in the white-noise cases. The points in Fig. [3] 
were obtained by a double limiting procedure. For a fixed 
resampling time ires, we increased the number of samples 
until convergence was reached. Then we iterated the 
scheme for increasing values of ^ros until a stable value 
for A2 was obtained. The larger the resampling time, the 
larger the number of samples to keep the error within the 
chosen bounds. 



VI. POISSON SHOT NOISE 

In a dilute gas with short-range interactions, phase- 
space coordinates evolve trivially in-between collisions. 
During collisions, positions remain essentially unchanged 
while momenta experience sudden jumps. The same de- 
scription applies to tangent-space coordinates. Thus, in 
a mean-field setting, the tangent dynamics of a repre- 
sentative (effective) particle is described by Eq. ([8)), the 
stochastic frequency corresponding to Poisson shot noise: 



K(i) = ^ (5(t - i,) . 



(33) 



Neglecting correlations among collisions the amplitudes 
Ai will be modeled by independent stochastic variables 
(identically distributed). Accordingly, the succession of 
collision times {ti} constitutes a Poisson process. 

The random oscillator (|S]) with Poisson frequency ([55)1 
was solved by van Kampen [s^l (including damping and 
additive noise). He derived an exact integro-differential 
equation for the probability distribution P{q,p,t) from 
where the evolution of the moments (q^p™) can be sis- 
tematically obtained [s^. For the second moments one 
gets 



2 \ 

p{A^) -2p{A) 

-p{A) 1 ; 




(34) 



where p is the collision frequency. 

Remarkably the expression above can be shown to co- 
incide with the result of the second-order cumulant ap- 
proach (|13ll4p . However, the higher-order cumulants of 
njt ) are not null, rather, they are delta-correlated [5Q|- 
I52I I . Just they do not affect the asymptotic growth of the 
second moments. 

The equations for the fourth moments (g'*), (p'*), 
{q^p"^), {q^p), {qp^) can also be calculated without much 
effort. The corresponding matrix reads 



/ 

p{A^) 6p(A2) 

p{A^) 

-p{A) 3 

-p{A^) 1 ~3p{A} 



4 \ 

-4p(A3) -4p{A) * 

-2p{A) 2 



3p(A2) / 



(35) 
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from where one extracts the fourth-order generahzed ex- 
ponent A4. 

Figure |4] shows that, when importance-samphng is 
used, the agreement between theory and numerics is ex- 
cehent. On the other side, simple samphng (plus a Gaus- 
sian approximation) leads to deviations from the theory, 
which become stronger as collision frequency ( "density" ) 
is lowered. Of course this disagreement is a consequence 
of the nonGaussianity of the distribution of finite-time 
Lyapunov exponents, and can also be observed when 
comparing A vs 2X2 — A4 (this can be thought of as a 
failure of the replica trick [l6|, [HI in its crudest version) . 
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time and number of samples) that result in a faster con- 
vergence. 

On the theoretical side, we carried out several tests of 
the cumulant approach in nontrivial cases, i.e., for fre- 
quencies corresponding to Ornstein-Uhlenbeck and Pois- 
son processes. In particular, we identifyed the cor- 
rect perturbative parameter (Kubo number) and -not 
unsurprisingly- verified that the second-order truncation 
of the cumulant series gives the exact second-order gen- 
eralized exponent A2 for the case of Poisson shot noise. 

Concerning the application of the cumulant approach 
to dilute gases, we note that in this case the tangent dy- 
namics can be thought to be driven by multivariate Pois- 
son noise. Accordingly the second-order truncation could 
indeed produce the exact A2 -like in the one-dimensional 
problem. However, the verification of this expectation 
would require the numerical calculation of A2 for a Hamil- 
tonian, i.e., nonstochastic, system. In order to imple- 
ment an importance-sampling algorithm for this case one 
should somehow introduce noise in the dynamics, then 
calculate A2 as a function of the noise intensity, and ex- 
trapolate the results to zero noise [s^. Several ideas for 
constructing such an algorithm are currently under in- 
vestigation. 

Acknowledgements: 



FIG. 4: (Color online.) Harmonic oscillator with Poisson- 
shot-noise frequency. We show the Lyapunov exponents A, 
A2 and A4 as a function of collision frequency p. Red/blue 
lines indicate theoretical estimates for A4/A2, and the corre- 
sponding symbols stand for numerical calculations using ei- 
ther simple-Gaussian sampling (open symbols) or importance- 
sampling (full symbols). Shown is also the theoretical result 
for 2A2 — A4 (black line) , which is an estimate of the standard 
Lyapunov exponent A (circles, numerical). 

The numerical method worked satisfactorily, the rela- 
tion between parameter values and efficiency being sim- 
ilar to the white- noise cases analyzed in Sect. IIVI 



VII. FINAL REMARKS 

We analysed the random harmonic oscillator as a sim- 
plified model of the tangent dynamics of many-particle 
systems. In spite of its relative simplicity, this model 
already exhibits some of the essential features and char- 
acteristics of high-dimensional systems. 

Specifically, we were able to assess the performance of 
the importance-sampling approach for the numerical cal- 
culation of generalized Lyapunov exponents. In all the 
considered cases -some of which unaccessible by stan- 
dard sampling methods- we confirmed that the method 
works satisfactorily, and developed some intuition about 
the appropriate values of the parameters (i.e., resampling 
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Appendix A: Fourth cumulant 

Here we briefly describe the calculation of the fourth- 
cumulant contribution to the generalized Lyapunov expo- 
nent of the Ornstein-Uhlenbeck oscillator, i.e., the right- 
most term in Eg . . In general, this contribution 
reads: [Illli,!!!! 

Ki{t) = K^-*) - K(2) = e^"'Q4t) e-^<>' , (Al) 

where 

rt nil rii 

Qi{t) ^ dti dt2 / dh 
Jo Jo Jo 

- (^Bi{t)Bi{h)) (^Bi{t2)Bi{h)) 

- (^Biit)Bi{t2)) (^Bi{ti)Bi{t3)) 

- (i?i(t)i?i(i3)) (Mti)Bi{t2))) , (A2) 

with 

Bi{t) ^ e-^"'Bi{t)e^'>* . (A3) 
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For the large times we are interested in, i.e., t ^ t, K/^it) 
becomes time independent. Next we note that matri- 
ces Bi are proportional to the scalar Ornstein-Ulhenbeck 
process £,{t). So, one must only calculate two- and four- 
time correlators of S,{t). By virtue of the Gaussian prop- 
erty [IB] , the four-time correlator is expressible as a sum 
of products of two-time functions ([28| . Finally one calcu- 
lates the triple integrals and takes the limit t ^ oo (with 
the help of an appropriate software, e.g., Mathematica 
[stI]). arriving thus at the desired result (pO|) . 



Appendix B: Lyapunov exponent 

Here we sketch the steps leading to the approximate 
expression for the Lyapunov exponent (of the random- 
frequency Ornstein-Ulhenbeck oscillator) that is plotted 
in Fig. [3l We have simply adapted the calculations of 
Mallick and Peyneau [40] to the case kq = 0. 

In the absence of damping, the Lyapunov exponent can 
be obtained as '40] 

A = lim l(ln(g2^a> (Bl) 

I— ^oo Zt 

= hm \±an{q^ + e)) (B2) 
t-!-oo 2 at 

= lim((2/) + i-|(hV + l))) (B3) 
t->co 2 at 

= lini(j/), (B4) 

where y = q/ q. From Eq. ^ one sees that y obeys the 
following nonlinear equation 

y^-y'' +T^{t). (B5) 



We will first find the exact expression for A when the 
noise is white (intensity A). In this case, the associated 
Fokker-Planck equation for P{y,t), i.e., 

dtP = dy{y'P) + ^dyyP, (B6) 

has the following steady state solution: 

P.M = Ne-^y"'^^'''^ r e^^'^^'^^dx, (B7) 

J —oo 

where is a normalization constant. By averaging over 
the steady state we obtain 

A-(A) = dyyPUv) = ^ (^) ' 0.2893 . 

(B8) 

In the case of an arbitrary correlation time r, by using 
a mean-field approximation ("decoupling ansatz" p^). 
one can derive the following equation for A: 

So, the final result comes in the form of an implicit equa- 
tion: 

A(A, r). 0.289 (^^^A__y. (BIO) 

This approximate relation slightly underestimates the 
numerical results of Fig. [3l 
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